function [ model,Theta2 ] = nntrain(X,y,hidden_layer_size,num_labels,lambda)
%    compute (input->hidden) theta1 and (hidden->output)theta2
%   [ Theta1,Theta2 ] = nntrain(X,y,hidden_layer_size,lambda)
%

input_layer_size  = size(X,2);  
if nargin<4
    num_labels = length(unique(y));         
end

m = size(X, 1);

if nargin<5
    lambda = 0;
end

initial_Theta1 = rands(hidden_layer_size,input_layer_size+1);
initial_Theta2 = rands(num_labels,hidden_layer_size+1);
initial_nn_params = [initial_Theta1(:) ; initial_Theta2(:)];



                   

options = optimset('MaxIter', 50);

costFunction = @(p) nnCostFunction(p, ...
                                   input_layer_size, ...
                                   hidden_layer_size, ...
                                   num_labels, X, y, lambda);
                                   
                                   
[nn_params, ~] = fmincg(costFunction, initial_nn_params, options);


Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), ...
                 hidden_layer_size, (input_layer_size + 1));

Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), ...
                 num_labels, (hidden_layer_size + 1));




model.theta1 = Theta1;
model.theta2 = Theta2;
model.lambda = lambda;

if nargout>1
    model = Theta1;
elseif nargout==1
    model.theta1 = Theta1;
    model.theta2 = Theta2;
    model.lambda = lambda;
end

end